Continuous time multidimensional Markovian description of Levy walks 
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The paper presents a multidimensional model for nonlinear Markovian random walks that gener- 
alizes one we developed previously (Phys. Rev. E 79, 011110, 2009) in order to describe the Levy 
type stochastic processes in terms of continuous trajectories of walker motion. This approach may 
open a way to treat Levy flights or Levy random walks in inhomogeneous media or systems with 
boundaries in the future. The proposed model assumes the velocity of a wandering particle to be 
affected by a linear friction and a nonlinear Langevin force whose intensity is proportional to the 
magnitude of the velocity for its large values. Based on the singular perturbation technique the 
corresponding Fokker-Planck equation is analyzed and the relationship between the system param- 
eters and the Levy exponent is found. Following actually the previous paper we demonstrate also 
that anomalously long displacements of the wandering particle are caused by extremely large fluc- 
tuations in the particle velocity whose duration is determined by the system parameters rather than 
the duration of the observation interval. In this way we overcome the problem of ascribing to Levy 
random walks non-Markov properties. 

PACS numbers: 05.40.Fb, 02.50.Ga, 02.50.Ey, 05.10.Gg, 05.40.+a 



I. INTRODUCTION 

There is a wide class of physical systems, where trans- 
port phenomena exhibit anomalous behavior of the Levy 
type called also superdiffusion (for a review, see, e.g., 
Ref. [l|). For such processes a variable x describing, in 
particular, spatial displacement of a wandering particle 
in the iV-dimensional space 1* during time interval t is 
characterized by some scale £(t) increasing with time t 
as £(t) oc t x l a , where the parameter a meets the inequal- 
ity a < 2. In addition, the distribution function P(x, t) 
of the particle displacement has the asymptotic behavior 
ln[P(x,<)] ~ -e{t) a /\x\ N+a for |x| »£(t). 

In particular, it is the case for the motion of tracer 
particles in turbulent flows [2| , the diffusion of particles 
in random mediaH , human travel behavior and spread- 
ing of epidemics [j] or economic time series in finance Q . 
Recently, there has been a great deal of research about 
superdiffusion. It includes, in particular, a rather general 
analysis of the Langevin equation with Levy noise (see, 
e.g., Ref. [6|) and the form of the corresponding Fokker- 
Planck equations @, H, U EH) description of anoma- 
lous diffusion with power law distributions of spatial and 
temporal steps fill . Levy flights in heterogeneous 
media [H, 0, Il5| and in external fields [lfl |T3|, con- 
structing the Fokker-Planck equation for Levy type pro- 
cesses in nonhomogeneous media 0, [H, Hof, first pas- 
sage time analysis and escaping problem for Levy flights 
HE HI HI IE HE HI HI , as well as processing exper- 
imental data for detecting the Levy type behavior [23 |. 
Besides, it should be noted that the attempt to consider 
Levy flights in bounded systems (see, e.g., Ref. [30l l3l| 
and review [l[ as well) has introduced the notion of Levy 
walks being a non-Markovian process because of the ne- 



cessity to bound the walker velocity. The problems arise 
because for a Levy walk the second moment ((x t — x^) 2 ) 
(xj and x t are the initial and final point, respectively) 
diverges. 

Previously, we developed an one-dimensional model 
that generates continuous Markovian trajectories, follow- 
ing the Levy statistics, by using simple Gaussian but 
multiplicative noise for the time evolution of the velocity 
[32| . The spatial dynamics naturally follows from this. 
It should be pointed out that a first step in this direction 
can be found in Refs. [HI, 34] . For a fixed time scale St we 
can recover the standard behavior of Levy type processes. 
However, we have full locality in the sense that a trajec- 
tory can be determined with any desired resolution. In 
other words, the developed model proposes a microscopic 
implementation of the Levy type processes characterized 
by an arbitrary small time scale r that can be chosen 
beforehand. When running time exceeds essentially this 
microscopic time scale, Or, the corresponding random 
walks, as should be, are described by the Levy distribu- 
tion. 

The purpose of the present paper is, first, to general- 
ize the developed one-dimensional model to multidimen- 
sional case and, second, to obtain the rigorous results 
using an original singular perturbation technique. The 
previous paper [32j was actually devoted to the formu- 
lation of the problem at hand and the qualitative expla- 
nation of its main properties using numerical simulation. 
In the present paper the one-dimensional model is also 
analyzed in detail as a specific limit case and the cor- 
responding results just declared previously arc obtained 
in a rigorous way. The approach to be developed may 
open in the future a way to treat Levy flights or Levy 
random walks in inhomogeneous media or systems with 
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boundaries including boundaries of complex geometry. 



Essence of the proposed description 



II. STOCHASTIC SYSTEM AND THE 
GOVERNING EQUATIONS 

We consider continuous random walks in the Eu- 
clidean A-dimensional space 1* governed by the fol- 
lowing stochastic differential equations of the Hanggi- 
Klimontovich type (a post-point process in the given 

case) mm 1| 



dxi 
~dt 

dvi 



(1) 



___{N±^l Vi + ^ g{v)Hi{t)) (2) 

for i = 1, 2, .... N. Here the symbol * denotes the be- 
longing of the given product to the Hanggi-Klimontovich 
type, x = {xiX2, • ■ ■ , xn} is the point of M. N , the vector 
v = {i>i,i>2) ■ ■ ■ j vn} is the current particle velocity, and 
the collection of mutually independent random variables 
feO)}i=f such that 

(&(*)> = and (£i(t)&(t'))=du>8(t-t') (3) 

represents the Langevin forces with the amplitude 

5 (v) = v^f+v^. (4) 

The dimensionlcss coefficient a > 0, the time scale r, and 
the characteristic velocity v a are the system parameters. 
According to the results to be obtained, the Levy type 
behavior can be found for 



1 < a < 2. 



(5) 



Exactly this region of a will be considered below. The 
scalar form of the coefficient a and the Langevin force 
intensity (j4|) actually implement the adopted assumption 
about the isotropy of the system at hand. For N = 1 the 
present model coincides with one developed in Ref. [32| 
within the replacement r — > 2t. 

For the given system the distribution function 
V (x — xo, v, vo, t) obeys the following Fokker-Planck 
equation written in the kinetic form 



dV 
dt 



N 

£ 

i=l 



f - — 

T dVi 



dV 



dvi 



subject to the initial condition at t = 

V(x- x , v, v, 0) = 5(r - r )5(v - v ) 



(6) 



(7) 



where, an addition, the system translation invariance 
with respect to the variable x is taken into account ex- 
plicitely. 



Let us discuss the core idea of the mathematical de- 
scription of the Levy type random walks within the con- 
tinuous time Markovian process. In the general form the 
corresponding Fokker-Planck equation ([6]) can be rewrit- 
ten as 



dP 
dt 



-v- V, 



V 



(8) 



where the operator L v acting on the variable v only is 
given by the second term on the right hand side of ex- 
pression ©. Since our analysis is confined to the long 
time dynamics of the variables x the solution V of this 
equation is sought in the form of the expansion 



P = ^>a(v,V x )/ a (x,;) 



(9) 



over the eigenfunctions \P\ (v, V x ) matching the suitably 
defined eigenvectors A (V x ) of the linear eigenvalue prob- 
lem 



A(V x )(v,V, 



-v-V x + L v # A (v,V 



(10) 



where the operator V x is treated as some formal param- 
eter. We expect that the long time dynamics will be 
described only by a few terms of such expansion. These 
eigenfunctions are assumed to be complete and form a 
basis, which is justified as will be seen below. Further- 
more, according to the general properties of such ran- 
dom walks the zeroth approximation of the eigenvalue 
problem (|10p determined by the replacement V x — > 
describes the velocity distribution. In this case the cor- 
responding quantities {A := A (V x — ► 0)} are eigenvalues 
of the Fokker-Planck operator L v . In the collection {A} 
there is one zero eigenvalue denoted by Ao = and the 
related eigenfunction is the stationary distribution of the 
velocity v. All the other eigenvalues A have negative real 
parts, so that the corresponding eigenvectors describe de- 
caying deviations from the stationary distribution. 

Inserting ansatz ([9]) into the Fokker-Planck equa- 
tion ||5J) and separating terms with different eigenfunc- 
tions we obtain the governing equations for the coeffi- 
cients /a(x, t) 



a/A(x,f) 

dt 



= A(V x )/ A (x,i). 



(11) 



The quantities / A (x, t) except for one matching the eigen- 
value Ao(V x ) originating from Ao = decay in time 
rather fast due to the fact that the corresponding eigen- 
values A < of the zeroth approximation are negative 
and split from zero by a finite gap. So the long time 
dynamics of the random variable x is governed by a gen- 
eralized diffusion equation 



d/o(x,t) 
dt 



A (V x )/ (x,i) 



(12) 
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It is the standard diffusion equation if the eigenvalue 
A (V X ) is an analytic function of the operator V x . When 
the functional Ao(V x ) is non-analytic equation (fT2|) de- 
scribes anomalous diffusion. 

The purpose of the present work is to demonstrate that 
the non-analytic behavior of the eigenvalue Ao(V x ) be- 
comes possible when the second moment of the stationary 
velocity distribution corresponding to the Fokker-Planck 
operator L v diverges. 



III. VELOCITY DISTRIBUTION 

It is the statistical properties of the particle velocity v 
that are responsible for the Levy type dynamics of the 
given particle. So the present section is devoted to their 
individual analysis. It should qualitatively explain the 
results to be obtained below and elucidate the mechanism 
via which the Levy type process arises. 

The distribution of particle velocity is specified by the 
partial distribution function 



P„(v,vo,i) 



dxP(x-x ,v,v ,i) (13) 



obeying the reduced Fokker-Planck equation 

JV 



dP v _ d 
^ Or, 



dt 



8P 

g 2 (v)^ + (N + a)viP v 
ov. 



with the initial condition 



P„(v,v,0)=<5(v-v ). 



(14) 



(15) 



These expressions stem directly from equation © and 
the initial condition ([7]) after integrating them over R w . 

Below in this section the velocity extremum distribu- 
tion will be analyzed based on the first passage time prob- 
lem. For this purpose we need another equation govern- 
ing the velocity distribution P v (v, vo, t) and acting on the 
initial velocity vo. Namely, it is the following backward 
Fokker-Planck equation of the Ito type 



dP v 
dt 



N 

E 

i=l 



5 2 (v 



o 2 p v 
1 dv L 



(N + a- 2)v 0i 



dP v 



(16) 



conjugated with the forward Fokker-Planck equation (Q3 
of the Hanggi-Klimontovich type and subjected to the 
same initial condition (fT5"|) (see, e.g., Ref. [38j|). 



Stationary distribution 

Equation (fl"4)) admits the stationary solution P®*(v) 
obeying the condition of zero value probability flux 



apst 

g 2 M^- + (N + a)v t P^ 







(17) 



and the normalization to unity 



/ dvPf(v) = l. (18) 
Equation (fT7|) and condition (|18p directly yield 



^ st (v) 



T[(N + a)/2] < 



Tr^r (a/2) [g(v)} N + a ' 



(19) 



where T(. . .) is the Gamma function. It should be noted 
that expression (jTTJ) reflects the fact that the given ran- 
dom process admits the detailed balance with respect to 
the particle velocity treated individually. For the expo- 
nent a belonging to interval ([5]) the first moment of the 
velocity magnitude v = |v| converges whereas the second 
one diverges. 

The found expression (|19p for the velocity distribution 
actually indicates that the particle displacement 5r ~ vSt 
in the space M. N during some time interval St is likely to 
exhibit also similar power-law distribution typical for the 
Levy type random processes. 



Moment dynamics 

The Levy flights and, partly, the Levy random walks 
are characterized by mutually independent succeeding 
steps in the particle displacement. This subsection an- 
alyzing the dynamics of velocity moments illustrates us 
that the proposed model does exhibit the loss of correla- 
tion in the particle velocity on time scales exceeding the 
value r. Thereby the partition of a particle trajectory 
into segments of duration dt ;g> r really can be regarded 
as a sequence of independent particle jumps. 

In order to analyze the time dependence of the veloc- 
ity moments, equation (| 14[) is multiplied, at first, with a 
general function w(v) and integrated over all the possi- 
ble values of the particle velocity. In this way we get the 
equality 



dt V 



N «2 
O LU 



N 



duj 



dvi 



where the symbol (. . .) means the standard averaging 
procedure, namely, 



(...) = J rfv(...)P„(v,v ,i). 

Then setting w(v) = Uj, |v|, ViVj (i ^ j), and v 2 the ob- 
tained equation is converted into the governing equations 
for the corresponding velocity moments 



r^ = -(N + a-2) <«,) 

d(W\ 



dt 



= -(a-l)(\v\) + H N (t) 



(20) 
(21) 
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where 



subject to the boundary condition 



K N (t) 



'2v 2 a P v (0,t) for AT = 1, 

(JV-lWJ-\ foriV>l, 



t d{viV 3 ) 
2 dt 

T d(v 2 ) 

2 



(JV + a - 2) , for i^j (22) 



= (2-a)(v 2 ) + 7V w 2 



(23) 



According to equations (f20|) and (|22|) . on time scales 
about r the particle forgets the direction of its initial 
motion and the motion along different axes becomes in- 
dependent. Equation (|2"Tj) demonstrates the fact that the 
first moment of the velocity v remains bounded during 
the system motion and attains its equilibrium value at 
the same time scales. This feature which will be used fur- 
ther in analyzing the spectral properties of equation © . 
Finally, equation (f2"3"| shows us that the system has to get 
its stationary state on temporal scales actually exceeding 
the parameter r independently of the initial velocity. 



Velocity maximum distribution 

Previously [32| we claimed and justified numerically 
the fact that the Levy type behavior of the analyzed ran- 
dom walks in the ID case on time scales t 3> r is caused 
by the properties of extreme fluctuations in the parti- 
cle velocity. Namely, an anomalously long displacement 
of the particle during time interval t is approximately 
determined by its motion during the spike in the veloc- 
ity fluctuations with the maximal amplitude. We will 
substantiate this statement also in the general case com- 
paring the results of the present subsection devoted to 
maximum statistics of the given random walks and the 
distribution function of particle spatial displacements to 
obtained further. 

Let us consider the first passage time problem for the 
model at hand. The probability -F(vo,i?,t) for the par- 
ticle with initial velocity Vo such that |vo| < i) to reach 
the sphere |v| = i3 in the velocity space for the first time 
at the instant t is directly described by the backward 
Fokker-Planck equation (fl6|) . In particular its Laplace 
transform 

/>oo 

F L (v ,ti,s) = / dte- st F(v Q ,tf,t) 
Jo 

obeys the equation (see, e.g., [HI]) 



tsFt, 



N 

£ 

»=i 



2, ,d 2 F L 

uu 0a 



(N 



dF L 

dv -, 



Fi(v ,1?,«)| vo=tf = 1. (25) 
Due to the symmetry of problem (|24|1 . (|25l) its solution 
-Fl(vo,i9, s) is a symmetrical function of the argument 
Vo and can be treated as a function Fl(vo,'&, s) of the 
magnitude 



vo 



In this terms equation (|2"4")) can be rewritten as 




tsF l =g 2 {v ) 
- (N ■ 



8 2 F L (N-l)dF L 



dv 2 



2)v 



dF L 
dv 



(26) 



where the function g 2 (vo) is given by formula (|H). 

The introduced first passage time probability is nec- 
essary to analyze the velocity maximum distribution. 
Namely we need the probability $(vo,i?, t) for the ve- 
locity pattern v(i) originating from the point Vo such 
that Vo < $ to the get the maximum |v| = *d during the 
time interval t is related to the probability F(vq, t) by 
the following expression |39[ 



*(v ,7?,t) = ~JdeF{v ,^H) 



(27) 



or for the Laplace transforms 



$ L (v ,tf,s) 



ld_ 



F L (v ,#,s) 



(28) 



(24) 



So the desired probability $(vo, i?, t) = &(vo,i), t) as well 
as its Laplace transform actually depends on the magni- 
tude vo of the vector vo rather then on its components 
individually. 

To examine the characteristic properties of the first 
passage time statistics let us consider the limit case 
st < 1 and d v a . Under these conditions in agree- 
ment with the results to be obtained we can regard 
the function Fi J (vo,'d, s) to be approximately constant 
Fl{vo^is) — F n (-d 7 s) inside some neighborhood Qo OI 
the origin, vq = 0, whose thickness is much larger than 
v a - In particular, for s — > it is the sphere |v| < itself 
and Foifi, s) = 1 by virtue of (j25|) . Within the neighbor- 
hood Qo equation (|2"6"|) can be integrated directly with 
respect to the formal variable /(i>o) := dF^/dvo using 
the standard parameter-variation method. In this way 
taking into account that /(0) = due to the system 
symmetry we obtain the expression 



and for w ^> w, 



dF L (v ,$,s) 
dv 



™fl, ( «,.). r (« r »)< 



2r ( 



AM 



(30) 



Expression (|30|) . first, enables us to estimate the size 
$_l(s) of the domain Qo- In fact, inside the domain Qo 
the inequality 



Fo(iM) > wo 



dF L (vo,0, s) 
dv 



=> TS-V- < 1 



has to hold, which allows us to set 



^( s ) ~ ( — I v <> > v o 



or converting to the time dependence 



(31) 



V a > «0 ■ 



(32) 



So the characteristic velocity scale characterizing the first 
passage time probability and aggregating its time depen- 
dence is i?(t). 

Second, in the case d ^> #l(s) there is a spherical 
layer » a « % « d(t) (for N = 1 it is a couple of 
domains) wherein the equality F^wo, s) ~ -Pb($> s ) 
holds whereas the derivative dF^/dvo scales with w as 
dF^/dvo oc Vq^ 1 . This asymptotic behavior can be ob- 
tained also analyzing the solution of equation (|24|) for 
\vq\ 3> w a where <7 2 (wo) — Wq. In this case equation 
admits two solutions of the form 



with 



F L (w ,i9,s) oc v, 



5i ~ a and g 2 



91,2 



rs 
a 



(33) 



The second solution is relevant to the function 
Fl(vq, u\ s) only within the crossover from Fl(vq, i?, s) oc 
Wq to F^wo,??, s) pa Fo(-#, s) and even in this region, i.e. 
$(i) the derivative dF^/Ovo is determined by its 



< 



asymptotics Fl(vo,$, s) oc Wq. For larger values of wo, 
i.e., wo ^> $(t) the first passage time distribution is given 
by the expression 



F L (vo,$,s)^(^y 



(34) 



taking into account the boundary condition (|25l) . So we 
can write 



dF L (v o ,0,s) 
dv Q 



a- 



(35) 



also for wo < 

Expressions (|30"f and (|35|) describe the same asymp- 
totic behavior of the function Fi(wo,u\ s). Thereby we 
can "glue" them together, obtaining the expression 



2aT(^) 1 (Va 



r(£)r(f) rs 



(36) 



that holds in the limit d 3> ^l(s). It should be noted that 
this procedure is the kernel of the singular perturbation 
technique which will be also used below. 

Expression (136|) immediately gives us the desired for- 
mula for the maximum distribution ^(wo, w\ s). Namely, 
by virtue of ([28]) . for wq < $l(s) and $ $l(s) we have 



$ L (w ,u\s) 



2a 2 r(^±a) 1 wg 

r(f)r(f) rs 2 



(37) 



Then restoring the time dependence of the extremum dis- 
tribution from its Laplace transform the asymptotic be- 
havior for $ 3> $(t) we get 



r(f)r(f)r^ 



(38) 



Finalizing the present subsection it is possible to draw 
the conclusion that for wo <C $(t) the maximum velocity 
distribution is described by a certain function 



with the asymptotics 
$o (0 = 



2a 2 T 



r(f)r(f) e +1 



(39) 



(40) 



Here the velocity scale d{t) is given by expression (|32j) . 
We remind that distribution (|39p describes the magni- 
tude of the velocity maximum attained during time inter- 
val t. For wq -C ${t) all the directions of particle motion 
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are equivalent. So for a multidimensional space, N > 1, 
the probability density of the maximum velocity attained 
during time interval t being equal to i9 is given by the 
function 



$(«b,0,t) = ^-$(UD,|0|,t) 



(41) 



where Sn is the area of the sphere of radius in M. N 



Sn 



r(f) 



In ID space this sphere is degenerated into two points 
±-$, so the corresponding probability is 



$(«o,±iM) = ~*(ix>,iM) 



(42) 



because of the symmetry of the velocity fluctuations. 



IV. GENERATING FUNCTION 

General relations 

The present section formally substantiates that 
model ([1]), ([2]) does exhibit the Levy type behavior on 
time scales t ^> r. For this purpose let us analyze the 
properties of the generating function introduced as fol- 
lows 



G(x, k, t) := ( exp <^ [(x - x ) x + rvk] 



(43) 



where the used averaging operator (...) is specified by 
the expression 

(...)= J djcdv(...)7 , (x-xo,v,v ,t) (44) 



and the wave vectors x and k are dimensionless variables. 
Then by virtue of the Fokker-Planck equation $Q and 
the initial condition ([7]) the generating function obeys 
the governing equation 



dt ^ \ dh \ dh 

%— 1 v x 



+ [xt- (N + o)h]^- \ -k 2 G (45) 



subject to the initial condition at t = 



Q(x, k, 0) = exp <^ — v k 



(46) 



At the origin k = and x — the function (|4"3"|) meets 
also the condition 



Q(o,o,t) = i 



(47) 



which follows directly from the normalization of the dis- 
tribution function to unity. In deriving equation (|4"5|) the 
following correspondence between the operators acting in 
the spaces {x, v} and {x, k} 



d 

dxi 



d_ 

dv. t 



-w, 



d_ 

a dh 



as well as the commutation rule 

Oh; 1 Kj dki~ ij 

have been used. 

The argument x enters equation (|45[) as a parameter; 
the given equation does not contain any differential oper- 
ator acting upon the function Q(x, k, t) via the argument 
x. This feature enables us to raise a question about the 
spectrum properties of equation (|4"S"|) , where the variable 
x is treated as a parameter given beforehand. Then the 
desired eigenfunctions and their eigenvalues 



{* A (k|^,p)}, {A(^,p)} 
obey the equation 



(48) 



A* A 



N 

E 



d_ 



dk t 



+ [xi -(N + a)h 



dh 



k 2 */ 



(49) 



where the symbol p denotes the complete collection of the 
eigenfunction parameters for a fixed value of x. We point 
out that the time dependence corresponding to these 
functions has been chosen in the form exp{— At/r}, which 
explains the minus sign on the left-hand side of (I45|) and 
dimensionless type of the eigenvalues {A (x, p)}. 

In these terms the solution of equation (|45[) can be 
written as the series 

G(x, k, i|v ) = ^f(x, p|v ) * A (k|x, p) 



Xexp|-A(*r,p)-j , (50) 

where {f(x, p|xo, vo)} are the coefficients of expan- 
sion (|50p which meet the equality 

f (x, p|v ) *a (k|x, p) = exp I — v k 1 (51) 



following from the initial condition (|46p . In agreement 
with the results to be obtained below the spectrum of the 
Fokker-Planck equation (|45|) is bounded from below by a 
nondegenerate minimal eigenvalue A mm (x) > 0, whereas 
the other eigenvalues are separated from it by a final 
gap. So, as time goes on and the inequality t ^> t holds, 
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Q{k, k, t) = f miD (x|v ) * min (k|x) exp ^ -A min (x 



the term corresponding to the minimal eigenvalue will be 
dominant and sum (1501) is reduced to 



t 



(52) 

on large time scales. Here ^min (k|x) is the eigenfunction 
matching the eigenvalue A mm . 

Whence several consequences follow. First, the general 
identity (|47p holds at any moment of time, thereby 



and due to J58 



A min (0) = . 



(53) 



Second, on large time scales f>r the system has to "for- 
get" the initial velocity Vo, so the expansion coefficient 
/min (x) does not depend on Vo and it can be aggregated 
into the function ty m i n (k|x). In this way the initial con- 
dition (I5T1) reads 



«W(k|x)+ J2 f ( x ^l v o) *A (k|x, p) 

p,A>A min 



exp <J — v k } (54) 



The terms in sum with A > A m j n determine the 
dependence of the generating function C/(x, k,t|vo) on 
the initial velocity vo and, thus, the corresponding coef- 
ficients / (x, p|vo) must depend on Vo. Finding the first 
derivative of both the sides of this equation with respect 
to v we have 



E 

P,A>Amin 



d/(x, p|v ) 
0v a 



*a (k|x, p) 



h exp \ — v k > (55) 



v a (. v a 

and for k = and any value of the initial velocity Vo 
<9/(x, p|v ) 



E 

3,A>A min 



dv ., 



-f A (0|x,p) = 0. (56) 



The latter feature, third, enables us to write individually 

* A (0|x,p) = 0, A>A min (57) 

for all the eigenfunctions except for "Jmin (k|x). By virtue 
of (|54"|) and condition (|57p the eigenfunction 'I'mjn (k|x) 
meets the equality 



*min (0|X) = 1 



(58) 



at k = 0. 

Summarizing the aforementioned we see that on large 
time scales t 3> r the asymptotic behavior of the given 
generating function is the following 

^(x, k, t) = * min (k|x) exp I -A min (x) - i , (59) 



G(x, 0, t) = exp \ -A min (x) - 



(60) 



for k = 0. In what follows calculating the eigenvalue 
A m in (x) will be the main goal. 

The randow walks under consideration should exhibit 
the Levy type behavior on large spatial and temporal 
scales, i.e. for |x — xo v a r and J » t. This allows 
us to confine our analysis to the limit of small values 
of x, i.e. to assume |x| <C 1 and also A m i n (x) 1. 
Under such conditions the spectrum problem (149p may 
be studies using perturbation technique with term 



JV 

E- 

i=l 



(61) 



playing the role of perturbation. 



Zeroth approximation. Spectral properties of the 
velocity distribution 



The zeroth approximation of equation (j49J) in pertur- 
bation (foTj) matches the case x = 0, where the generating 
function G(0, k, i) actually describes the velocity distri- 
bution as stems from its definition (|43p . Setting x = in 
the eigenvalue equation (14"9")) we reduce it to the following 



A' 



(N + a)h 



where the designations 

$ A (k|p) = ¥ A (k|0,p) and A(p) = A(0,p) (63) 

have been used. 

Pursuing different goals let us consider the conversion 
of equation (|62p under the replacement 

<f A (k|p) = |k|^ A ,„(k|p) (64) 
for several values of the exponent (3 n . 

Case 1: /3 1 = (N + a)/2 

In this case the substitution of into equation 
converts it into 



N 







a 2 - N 2 



5A,1 



(65) 
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The operator on the right-hand side of equation (|65|) is 
Hermitian within the standard definition of scalar prod- 
uct. So all the eigenvalues {A} are real numbers and 
the corresponding cigenfunctions {4>X,n (k|p)} form a ba- 
sis. It should be noted that the given conclusion coin- 
cides with the well known property of the Fokker-Planck 
equation for Markovian systems with the detailed bal- 
ance [40j . In addition the eigenfunctions (p\,i (k|p) can 
be constructed in such a way that the identity 



dk^l i (k|p) <j)y A (k|p') = 5 iz 



(66) 



hold for all of them except for the case describing the 
normalization of the eigenfunction 0o,i (k) corresponding 
to the minimal eigenvalue A m i n = A m i n (0) = by virtue 
of ([55]) . We note that the latter eigenfunction matches 
the stationary distribution (fT9")) of the particle velocity 
and its normalization is determined by condition (|58[) . 



Case 2: = (N + a - 2)/2 
In this case equation (|62[) can be rewritten as 



the value k = |k|. Indeed, expression (j68|) enables us to 
rewrite equation (|6"T|) as follows 



A0 



1 d 



A, 2 



N-l ' 



)\,2 



k N ~ 3 dk \" dk 
a 2 ~{N-2) 2 



k z + 



0a,2+L 2 a , 2 . (69) 



Thereby in the general case any eigenfunction <&a (k|p) 
can be written as the product of the corresponding sym- 
metrical eigenfunction and a certain eigenfunction of the 
angular momentum operator L 2 

= $ A (fc| Pa ).e w , 

and their eigenvalues are related as 

\(p s \Ju) = X(p s ) + SX u . 

So the given spectrum problem is splitted independently 
into its analysis with respect to the symmetrical eigen- 
functions and the spectrum problem for the angular mo- 
mentum operator. 



A0a,2 = -k 



N 

2 E 

i=l 



9A.2 



8B 



a 2 - (N-2f 



0A,2 ■ (67) 



Let us split the Laplacian entering the right-hand side 
of expression (|67p into two parts acting either on the 
magnitude k := |k| of the vector k or on its angular 
variables 



d 2 



N 

2—1 L 



i o 

k^dk 



kN 'ik) P L2 ' 



(68) 



where L is the angular momentum operator. We will re- 
mind the particular expressions of L 2 for two- and three- 
dimensional spaces 



T 2 

^=2 — 



T 2 



d 2 

d<p 2 ' 
1 d 



sin0- 



88 



sin 2 9 dip 2 



In one-dimensional space the angular momentum opera- 
tor li% = i can be treated as the symmetry operator under 
the reflection x \ — > —x. All the eigenfunctions {O w } of 

the operator L 2 possess eigenvalues {5X U > 8\$} exceed- 
ing some positive number <5Ao > of order unity, nat- 
urally, except for the function l w not depending on the 
angular variables for which L 2 1 W = 0. 

This feature enables us to confine our further analy- 
sis to the eigenfunctions {$ A (k\ps)} depending only on 



Case 3: = a/2 

For the given value of the exponent f3 and the symmet- 
rical eigenfunctions 4>\ 3 (k\p s ) equation (|6"2|) is reduced 
to the modified Bessel differential equation 



'A,3 



dej) 



A. 3 



dk 2 



dk 



a" 
T 



^A,3 



. (70) 



Since the desired eigenfunctions should decrease as k — > 
oo the solution of equation (|70|) is given by the modified 
Bessel function of the second kind 



4>U(k) ^ K v (k) 

with the order v = \J a 2 /4 — A because 



(71) 



K v {k) ~ *l — e as k — > oo . 
y 

for any value of the parameter v (4lj . 

Whence it follows that, first, there are no eigenfunc- 
tions with eigenvalues A < 0. Indeed, for small values 
of the argument k the modified Bessel function K v (k) oc 
k~ v when its order is a positive number v > 0, which is 
the case for A < a 2 /A. Under these conditions the trial 
function 

$(fc) := k^K v (k) cx k~ v+a/2 for k < 1 

exhibits strong divergence with k — > 0. Second, as it 
must, the value A m ; n = is the minimal eigenvalue and 
the given trial function takes some finite value at k = 
0. The corresponding eigenfunction will be analyzed in 
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detail below. Third, the interval < A < a 2 /4 does not 
contain any additional eigenvalue. In fact, otherwise, the 
trial function 

0A,i(fc) = Mk)k p *- Pl oc k- N ' 2 K v {k) 

would give rise to strong divergence in the normalization 
condition ([M]) . Finally, the eigenvalues A > a 2 /A form 
the continuous spectrum of the given problem. The cor- 
responding eigenfunctions via the modified Bessel func- 
tion of the second kind with pure imaginary order exhibit 
strongly oscillatory behavior as k — > and due to the 
preceding cofactor k^ 3 meet condition (|57j) . 

In order to construct the desired eigenfunction <& m i n (k) 
matching the minimal eigenvalue A m ; n = let us make 
use of the expansion of the function K v {k) for small val- 
ues of the argument k 



K v {k) 



1 



r(i 



i> + i) 



0(k 2 ) 



(72) 



which is justified for the order < v < 1 (see, e.g., 
Ref. dH). For A = the order v — a/2 and the lat- 
ter inequality holds due to the adopted assumption ([5|) 
about the parameter a. Expression (|64|) and the obtained 
result ([7T|) specify the dependence of the eigenfunction 
^ > min(fc) oc k a l 2 K a / 2 (k) on the argument k. In the case 
under consideration the general condition (I58p reading 
^min(O) = together with asymptotics (|T2"|) enables us 
to find the preceding constant. In this way the desired 
expression 



$min(k) 



2— 

r(f) 



1 



k\ a T(^: 



+ 0(k 2 ) . (73) 



is got. Expression (f73|) finalizes the analysis of the zeroth 
approximation. 

Summarizing the aforementioned we draw the conclu- 
sion that at x = the spectrum of the Fokker-Planck 
equation (|4"5|) for the generating function (|4*3")l does con- 
tain zero eigenvalue A m i n (0) = corresponding to eigen- 
function (|73[) which is separated from higher eigenvalues 
by a gap equal to a 2 /A (in units of t). We note that 
the given statement is in agreement with the conclusion 
about the spectrum properties for a similar stochastic 
process with multiplicative noise [13, EH, E3 • 



The eigenvalue A m i n (x) for |x| <C 1. Singular 
perturbation technique 

On time scales t ^> r only rather small values of the 
wave vector x contribute substantially to the distribution 
of the random variable x, Indeed, according to (|59| or 
(1601) only the wave vectors x that meet the inequality 

AaunO)- < 1 
r 



asymptotics 
of solution 
of equation (76) 







dki 



j asymptotics of 
(~] zeroth approximation, 
' eigenvalue equation (62) 



region of gluing 
the asymptotics 



FIG. 1: Illustration of the singular perturbation technique in 
gluing the asymptotics the eigenfunction *& m in(fc) constructed 
in different limit region of the velocity wave number k. 

are essential. Due to A min (x) — > as x — > (see 
and the estimate A m j n ~ 1 for |x| ~ 1 this inequality is 
converted into |x| < >fi <gc 1 with x t — > as t — > oo. So 
in what follows the inequality |x| <C 1 will be assumed 
to hold beforehand. 

The perturbation term (|61| mixes the eigenfunctions 
of zeroth approximation affecting the spectrum property 
of the generating function. However in the case under 
consideration perturbation (|6ip disturbs the eigenfunc- 
tion $ min (k) with A m ; n = substantially only in a cer- 
tain small neighborhood of the origin k = wherein 
k| ^5 \ K \- ^ stems directly from the form of the eigen- 
function equation (j49| . Outside the neighborhood 
the eigenvalue function 1 S , m in(k|x) should practically co- 
incide with $ m j n (k). 

So in the space of the wave numbers x there is a spher- 
ical layer 



< Ikl < 1 



(74) 



wherein formula (|73|) in the limit fc 1 approximates not 
only the eigenfunction $ m i n (k) but also the eigenfunction 
*min(k|>r) (Fig. [1]). In addition within this layer as 
well as the neighborhood Q M , first, the the two eigen- 
value functions are practically equal to unity, \& ro i n (k|x), 
^mii^k) ~ 1. Second, any power term whose exponent 
exceeds a < 2 (see inequality ([5])) may be ignored. The 
latter concerns also the last term on the right-hand side 
of equation (|4*9"j) . 

Keeping the aforesaid in mind and splitting the func- 
tion W m i n (k|>f) into two parts 



* mln (k|x) = l-^(k|x). 



(75) 



where ^£>(k|x) <C 1 the eigenvalue problem under consid- 
eration is reduced to solving the equation 



N 

E 



d / k2 cv 



dki 



(N + a)ki 



dtp 
dk~. 



subject to the requirement of the solution exhibiting the 
asymptotic behavior 



</>(k|x) 



2) r(^) 



(77) 



for |x| <C |k| <C 1. In some sense expression ((77)) is a 
"boundary condition" of equation (|76|) . 
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Scaling relations 



be considered separately. In this way we get for £ < 



Using the scaling transformations 

h = \x\ Q, A min = \x\ a A min , ^(k|x) ^ |x|XC|n) 

(78) 

the eigenvalue problem ([76jl , ([77f is reduced to the equa- 
tion 



A 



N 

mm = 2_j 

i=l 



d 



dip 



(N 



(79) 

where n = {rii} is the unit vector parallel to the wave 
vector k and the asymptotic behavior (|77[) is converted 
into the asymptotics 



^(C|n)-(f) Q ^| as 



(80) 



at infinity. In some sense condition (|50"|) "glues" the 
asymptotic behavior of the eigenfunction \P m in(k|x) re- 
sulting from its properties for sufficiently large values of 
k together with one stemming from small values of k, 
in this case, specified by the solution of equation (|79")l . 
Exactly such a procedure is the essence of the singular 
perturbation technique. 

The eigenvalue problem ([79"]) , (f8T)]) does not contain any 
external parameter, enabling us to expect that A mlVL is a 
value about unity, vl m m ~ 1 • This estimate actually is the 
basic result of the present paper because it immediately 
leads us to the conclusion about the scaling properties of 
the random walks in the space that are governed by 
the system of equations ([J)) and ([2]). 

However before discussing the obtained results let us 
find the specific values of the eigenvalue Amin, with 
one-dimensional and multi-dimensional spaces being an- 
alyzed individually. 



One- dimensional model 



For the one-dimensional space 1 
becomes 



(N = 1) equation J7HD 



(81) 



Here we have assumed the unit vector n to be positively 
directed along the C _ axis and omitted the index i at the 
variable £. 

Equation (|8Tj) can be solved directly with respect to the 
variable dip/d( using the standard parameter- variation 
method. However this solution can exhibit a strong sin- 
gularity at £ = 0, so the regions £ > and C < have to 



T C = lcl exp [c 



C- 



and for £ > 

— = c r e x P , 

d( VC 



-i/C 



At 



(82) 



a 



(83) 



where the constants C±oo specify the asymptotic behav- 
ior of the derivative 

Thereby according to condition (180)) that should be im- 
posed on the asymptotic behavior of the solution ip(£) 
and converting into the corresponding asymptotic behav- 
ior of its derivative 



C+c 



-c_ 



2 a T 



'2 + Q N 



(84) 



Solution ([83]) diverges as ( — > unless the equality 



a 



+00 



A 



, mil , r _i e- 4 <ie=o 





holds, whence taking into account (|84|) we find the desired 
expression for the eigenvalue A m i n 



r(^) 



2«-ir(§)r(a) 



(85) 



In obtaining ()85j) the recurrence formulas for the Gamma 
function have been used. 



Multi- dimensional model 

For the multi-dimensional space M (N > 1) the so- 
lution of equation ([75)) subject to the asymptotic behav- 
ior (|80p turns out to depend only on two variables, the 
absolute value £ := |£| of the wave vector £ and the an- 
gular variable 



E (l«l 

i=l s 



(86) 
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This statement follows directly from the symmetry of 
equation ([75]) and the "boundary condition" ([50)) . Us- 
ing the variables {£, 9} we can write 



_d_ 



C d( 



d_ 

09 



and. thus, 



N q 



c 



d 



N 



d_ 

d (l-6> 2 ) 9 



C 96> 



assuming the given operators to act on functions of the 
arguments C an d 9 only. Then equation ([75)) can be 
rewritten as 



d( 2 



X =C : -7-7 -(a-l)C^ 



+ 7(l-* 2 ) ^ 



9(9 



9C 



(87) 



Thereby in the given case the eigenvalue problem is re- 
duced to solving equation ([87)) subject to the boundary 
condition ([50)) with respect to the function ip(£,0). 

The boundary 9 = ±1 is artificial; it arises via the 
change of variables ([55)) . So the solution should be an 
analytical function of the variable 9 also all the points of 
the closed interval 9 G [—1,1]. The latter feature enables 
us to seek the solution of equation ([57)) as a power series 
with respect to 9 



Substituting ([88)) into ([57)) and gathering terms with the 
same power of 9 we obtain the collection of individual 
equations for the components of expansion ([55)) , namely, 



,2 d 2 ^o 



(a - 1) A = -2^ 2 - + A n 



dC, 2 v ~ d( ~ T * C 



and for n > 1 



C -JPT- - {a- l)C-77- - n{N + n- 2)V>„ 



d( 2 



d( 

in + 1) 

= -(n+ 2)(n + l)V'n+2 ^ — V'n+i 

-tpn-i ■ (89b) 



c 



Let us, first, obtain the particular expression for the 
eigenvalue A m i n in the limit case N 3> 1. According 



to equation (j89b[) at the leading order in the small pa- 
rameter 1/N 



i dMO 



(90) 



whereas the other components are of higher order in 1/N. 
So within the given accuracy equation (|89ap becomes. 



2 d 2 ipo 
d( 2 



1 

NX 



(a-l)C 



(91) 



Then imposing the "boundary" condition (I80|) on equa- 
tion (|91|) and using the method of variation of parameters 
it is solved, giving us the expression 



dMO , a -i r 1 

— Q exp 



V 2 7 -2 2 AT2,l min 



2N( 2 



2 «-ir(f) 



(92) 



Therefor the function ^o(C) does not exhibit a singularity 
as C -> if 



2 4r 



v 2 y 



3a 

2-r2(|) 



(93) 



which gives us the desired estimate of the eigenvalue yl m i n 
at the leading order in 1/N. 

For an arbitrary value of N the system ([59")) could be 
analyzed numerically using, for example, the following 
algorithm. As follows from equations ([59")) . in a small 
neighborhood of the origin £ = the functions {V'n(C)} 
can be expanded into power series 

M0 = C- [A n0 + A nl ( 2 + A n2 e . . .] , 

where the coefficients {A nm } can be found by solving the 
system (|89l) . In particular, beforehand, it is possible to 
set A 00 = due to the general equality ([58)) and defini- 
tion ([75)1 of the function tp((, 9). In addition, the equality 
Aio = A m i n stems directly from equation (|89a[) . There- 
fore the given system of equations can be subjected to 
the following "initial conditions" 



^n| C=0 = 
#n 

d( C =0 
dipt 

d( C =o 



for all the values of the index n, 
= for n ^ 1, 



A min for n = 1. 



In this way the system of equations (|89a[) and (j89b[) is 
converted into an initial value problem for any given value 
of the parameter A m i n . The solution is proportional to 
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the quantity yl m i n = 1 so, af first, setting, for example, 
-^min = 1 an d solving this system numerically one finds a 
coefficient ui* of the asymptotics V'o(C) ~ W *C" as C ~~ * 00 • 
Then 



-^min — 



w*2°r(^) 



is the desired eigenvalue. 



V. LEVY TYPE DISTRIBUTION OF SPATIAL 
STEPS 

Summarizing the results of the previous section we 
draw the conclusion that the continuous time model fTJ) , 
(|2|) for the random walks governed by the given multi- 
plicative noise does exhibit the Levy type behavior on 
time scales exceeding essentially the parameter r. 

To justify this statement in detail, we take into account 
the scaling relations (|78|) and rewrite the generating func- 
tion (|6"0|) for the particle displacements in the space {x} 
in the form 



Q(x,0,t) = exp<^ -A 



t 



(94) 



which holds for f>r. Introducing the distribution func- 
tion Pj(x — Xo,t) of the particle displacements by the 
integral 



P x (x-Xo,t) = / dvP(x-x , 



,v ,t) (95) 



the generating function (f9"4"|) can be represented as 
Q(x,0,t) = J dxP x (x-x ,t)exp|^— (x — x ) • x| . 

R N 

(96) 

Comparing expressions (|94|) and (|96[) we see that, first, 
the distribution function P x (x — xo, t) = P x (fe, i) is sym- 
metrical, i.e. it should depend only on the length of the 
displacement vector Sx := |x — xq|. 

Second, it does not depend on the initial velocity vo 
of the particle, so on such time scales the description of 
these random walks can be confined to the spatial vari- 
ables only. 

Third feature is the time dependence of the spatial 
scale £(t) characterizing the particle displacements dur- 
ing the time interval t. To find it we assume the argu- 
ment of the generating function (|94[) to meet the estimate 
\x\ a {t/r) ~ 1 for values of the wave vector x such that 
(£(t)|x|) / (v a T) ~ 1, see formula (|96|) . When it immedi- 
ately follows that 



£{t) = (v a r) 



(97) 



Due to the adopted inequality the exponent l/a > 
1/2, so the time dependence of the characteristic particle 



displacement £(t) really describes a Levy type stochastic 
process. 

Fourth characteristics is the asymptotic behavior ex- 
hibited by the distribution P x (5x,t) of the particle dis- 
placements. This asymptotics matches spatial scales 
8x 3> £(t) and, correspondingly, small values of the argu- 
ment of the generating function (f9"4")) . Under these con- 
ditions expression is reduced to the following 



/1„ 



t 



2 / dxP x (Sx,t) sin 2 

R N 



(X - Xp) • X 

2v a r 



(98) 



obtained by subtracting the expression combined 
and (|96| from itself calculated at x = 0. In agreement 
with the results to be found below, the distribution func- 
tion P x (Sx, t) on scales Sx 3> £(t) is to exhibit the follow- 
ing asymptotic behavior 



P x (Sx,t) ~C 



WW 

(Sx) N + a 



(99) 



where the constant C can be found using expression 
Namely, asymptotics (|99[) enables us to rewrite expres- 
sion (|9"5)) for the one-dimensional case (N = 1) as 



A r 



t 



4C[£(t)] a 



d(Sx) 
{Sx) 1 + a 



(Sx)\x\ 



2v a r 



and for the multidimensional case (N > 1) as 

7T OO 

,„ t 



(100a) 



(Sx 



(Sx) 



N+a 



2C[£(t)] a JdOj d{5x) 



Sn-i{Sx sin 0) sin 2 



(&e)|x| cos© 



2v a r 



(100b) 



where 



n — 

2ir 2 r 



N-l 



r(f) 



is the surface area of the ./V-dimensional sphere. Calcu- 
lating integrals (|100p gives us the required values of the 
coefficient C for N = 1 



C = 



,sin(^)r(^) 

2<*-%r(f) 



(101a) 



where in addition expression (|85[) has been taken into 
account and for N > 1 



C = 



-A min sin r (s±i) r (^) . (ioib) 



Formula (|9~?| and asymptotics (|M)l are the characteris- 
tic features of random walks belonging to the Levy type 
stochastic processes. 
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Besides, comparing asymptotics (|99|) of the distribu- 
tion function P x {5x,t) and asymptotics (jlTj) of the ve- 
locity extremum distribution ^(vq, t) we see that these 
dependencies coincide with each other within the replace- 
ment Sx = p&T, where the cofactor p ~ 1 is a certain 
number of order unity. The value r characterizes the 
time correlations in the velocity fluctuations of particle 
motion, see formulae (f2"0" )) - (j2"3")) . Therefore anomalously 
long spatial jumps of the particle on time scales t ^> r 
can be treated as the displacements of this particle gained 
within the spikes of the extremum velocity fluctuations 
during the time interval t. 

Governing equation of superdiffusion 

According to expression (f9"4"| the generating function 
G x (x, t) := Q (x, 0, t) obeys the following equation 

T ^r = - A *a» H Q G x . (102) 
at 

By virtue of this generating function is no more 
than the Fourier transform of the distribution func- 
tion P x {yi — xo,i) describing the spatial displacements 
of wandering particles during time interval t in units of 
v a r. So returning to the original spatial variables x equa- 
tion (|102[) is converted into the following 

9 -§ = -° (-VlT /2 P, ■ (103) 
Here the quantity 

<T = Anin^T"" 1 (104) 

can be regarded as the coefhcient of superdiffusion and 
the operator entering the right-hand side of this expres- 
sion is the fractional Laplacian, being an integral oper- 
ator with kernel proportional to the right hand side of 
expression (|99p . For the rigorous construction of the 
fractional Laplacian and the corresponding details in 
its possible representations a reader is referred, e.g., to 
Ref. [H, Hi| . Formula (|104|) is actually the desired im- 
plementation of the governing equation (|12|) for superdif- 
fusion. 



VI. CONCLUSION 

Via Eqs. (fTJ) , ||2J| we have presented a model, which im- 
plements Levy flights on a "microscopic" level. In par- 
ticular this allows one to describe the trajectory, char- 
acterized by Levy statistics, in a continuous fashion on 
every given time scale St by choosing t <C St. Indeed, 
fixing any small duration St of the Levy walker steps we 
can choose the time scale r of the given model such that 
St 3> r and, as a result, receive the Levy statistics for the 
corresponding spatial steps on time scale St. Of course, 
the Langevin equation has to be updated on a time scale 



St <C t. Moreover, expressions (|L03[) and (|104[) demon- 
strate the equivalence of systems of given a oc w"r Q_1 . 
with respect to their asymptotic behavior. Thus, all the 
details of the microscopic implementation of Levy flights 
are expressed by the exponent a and the superdiffusion 
coefficient a. 

In our previous work we have derived the Levy behav- 
ior based on the numerically derived observation that 
|Aa;| oc d which translates the extremal behavior of d to 
that of Sx. Furthermore, we have restricted us to the ID 
case. In contrast, in the present work we have strictly 
calculated the lowest eigenvalue of the corresponding 
Fokker-Planck equation and thus the generating function 
for the particle displacement for arbitrary dimension. As 
outlined in this work this calculation made use of the 
singular perturbation technique. In particular, it is pos- 
sible to prove this proportionality and, furthermore, to 
calculate the proportionality constant by comparing the 
distribution functions (j4Tj) (or (|42|) for the ID case) and 
with each other. Thus, this model has finally found 
a mathematically strict solution. 

In physical terms the non-diffusive behavior enters 
by the multiplicative noise, which gives rise to a self- 
acceleration of the system and is intrinsically connected 
to nonequilibrium situations. More generally, the emer- 
gence of slow-speed and fast-speed periods, as generated 
by the present model, is of current interest in optimal ran- 
dom search theory and in the analysis of animal move- 
ment patterns (see, e.g., [U 51]). In these cases the 
searching phases tend to be associated with slow speeds 
whilst relocation phases tend to be associated with high 
speeds. Furthermore, the present approach can be inter- 
preted as a generalization of the Kramers- Fokker-Planck 
equation describing diffusion of particles, where via the 
choice g = const the noise is purely additive. 

Let us reiterate the possible applications of the present 
model, in particular for more than one dimension. (I) 
From a numerical perspective it is possible to generate 
a Levy flight based on a straightforward simulation of 
the Langevin equations ©. In the well-defined limit of 
small St the trajectory can be constructed with arbitrary 
precision. (II) The present approach allows the consider- 
ation of Levy flights together with boundary conditions 
or for heterogeneous media, e.g., by introducing a de- 
pendence of v a to depend on the location of the Levy 
walker. This is possible due to the strict locality and the 
Markovian behavior of our model. (Ill) The developed 
approach also opens an way to constructing the path in- 
tegrals for the Levy random walks based in the Wiener 
measure and, then, to develop a description of nonlinear 
Levy processes in a self-consistent way. 
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